//	--------------------------------------------------------------------
//	Copyright(C) 2006,2007 Zhao Yukun. All rights reserved.
//	This file is a part of the Xeres system.
//	Xeres 2007
//	Contact : xeres.engine@gmail.com
//	$(reservedInfo)
//	$(reservedDeclaration)
//	--------------------------------------------------------------------
///	@file	<math_ext.h>
///	@path	~/src/lib/math/
///	@date	2007/09/08
///	@desc	mathematics extension.

#pragma once

// CRT
#include <limits.h>
#include <math.h>

// half
#include "lib/math/half.h"

namespace xeres {

	// for vector math
	enum AXIS
	{
		AXIS_X = 0,
		AXIS_Y = 1,
		AXIS_Z = 2,
		AXIS_last ,
	};

#define	CHECK_INDEX( INDEX , LIMIT )		( assert( (INDEX >= 0) && (INDEX<LIMIT) ) )

	// PI value
#define	MATH_PI				(3.141592653589793238462643383279502884197169399375105820974944592308)
#define	MATH_DOUBLE_PI		( 2.0 * MATH_PI )
#define	MATH_HALF_PI		( MATH_PI / 2.0 )
#define	MATH_PI_OVER_180	( MATH_PI / 180.0 )
#define	MATH_PI_UNDER_180	( 180.0 / MATH_PI )

	// legal arithmetic operand type trail
	// ( only int , float , double is legal ! )
	template< typename _Ty > struct ArithmeticOp;

	// int
	template<> struct ArithmeticOp< int >
	{
		// is math op
		static void verify( void ) {}
		// epsilon
		static int epsilon( void ) {
			return 0;
		}
		// comparison epsilon
		static int comparison_epsilon( void ) {
			return 0;
		}
		// no Pi functions
		// infinite
		static int infinite( void ) {
			return INT_MAX;
		}
		// negative infinite
		static int ninfinite( void ) {
			return -INT_MAX;
		}
	};
	// int
	template<> struct ArithmeticOp< uint >
	{
		// is math op
		static void verify( void ) {}
		// epsilon
		static uint epsilon( void ) {
			return 0;
		}
		// comparison epsilon
		static uint comparison_epsilon( void ) {
			return 0;
		}
		// no Pi functions
		// infinite
		static uint infinite( void ) {
			return UINT_MAX;
		}
		// negative infinite
		static uint ninfinite( void ) {
			return 0;
		}
	};
	// float
	template<> struct ArithmeticOp< float >
	{
		// is math op
		static void verify( void ) {}
		// epsilon
		static float epsilon( void ) {
			return FLT_EPSILON;
		}
		// comparison epsilon
		static float comparison_epsilon( void ) {
			return 1.0e-5f;
		}
		// PI
		static float pi( void ) {
			return float( MATH_PI );
		}
		// pi * 2 
		static float double_pi( void ) {
			return float( MATH_DOUBLE_PI );
		}
		// half pi
		static float half_pi( void ) {
			return float( MATH_HALF_PI );
		}
		// pi over 180
		static float pi_over_180( void ) {
			return float( MATH_PI_OVER_180 );
		}
		// pi under 180
		static float pi_under_180( void ) {
			return float( MATH_PI_UNDER_180 );
		}
		// infinite
		static float infinite( void ) {
			return FLT_MAX;
		}
		// negative infinite
		static float ninfinite( void ) {
			return -FLT_MAX;
		}
		// get NaN
		static float nan( void ) {
			static const int a = 0x7f800000;
			static const float * f = (float*)(&a);
			return *f;
		}
	};
	// double
	template<> struct ArithmeticOp< double >
	{
		// is math op
		static void verify( void ) {}
		// epsilon
		static double epsilon( void ) {
			return DBL_EPSILON;
		}
		// comparison epsilon
		static double comparison_epsilon( void ) {
			return 1.0e-12;
		}
		// PI
		static double pi( void ) {
			return MATH_PI;
		}
		// pi * 2 
		static double double_pi( void ) {
			return MATH_DOUBLE_PI;
		}
		// half pi
		static double half_pi( void ) {
			return MATH_HALF_PI;
		}
		// pi over 180
		static double pi_over_180( void ) {
			return MATH_PI_OVER_180;
		}
		// pi under 180
		static double pi_under_180( void ) {
			return MATH_PI_UNDER_180;
		}
		// infinite
		static double infinite( void ) {
			return DBL_MAX;
		}
		// negative infinite
		static double ninfinite( void ) {
			return -DBL_MAX;
		}
		// get NaN
		static double nan( void ) {
			static const __int64 a = 0x7ff0000000000000LU;
			static const double * f = (double*)(&a);
			return *f;
		}
	};

	/*
		\struct	MustBeFloat
		\brief	Verification of operation must be a float value.
	*/
	template< typename _Ty > struct MustBeFloat;

	template<> struct MustBeFloat< float > { static bool verify( void ) { return true; } };
	template<> struct MustBeFloat< double > { static bool verify( void ) { return true; } };
	template<> struct MustBeFloat< long double > { static bool verify( void ) { return true; } };

#define	TYPE_MUST_BE_FLOAT( TYPE )		assert( MustBeFloat<TYPE>::verify() );

	/*!
		proper type-select templates!
		float_cast - cast to floating point type ,
									e.g. int - float , double - double
		float_ar_select - float arithmetic operation select ,
									e.g. int & float will select float , float & double will select float
		precision_cmp - precision comparison , low precision will always at left
									e.g. float & int will be int at left , float at right
		precision_select - desire high precision
									e.g. float & double will select double
	!*/

	// float cast
	template< typename TT > struct float_cast
	{
		typedef float result;
	};
	// double
	template<> struct float_cast<double>
	{
		typedef double	result;
	};
	// long double
	template<> struct float_cast<long double>
	{
		typedef long double result;
	};

	// precision comparison
	template< typename T1 , typename T2 > struct precision_cmp;
	// the same type
	template< typename TT > struct precision_cmp< TT , TT > 
	{
		typedef TT left;
		typedef TT right;
	};
	// int & float
	struct precision_cmp_int_float
	{
		typedef int left;
		typedef float right;
	};
	template<> struct precision_cmp< int , float > : public precision_cmp_int_float {};
	template<> struct precision_cmp< float , int > : public precision_cmp_int_float {};
	// int & double
	struct precision_cmp_int_double
	{
		typedef int left;
		typedef double right;
	};
	template<> struct precision_cmp< int , double > : public precision_cmp_int_double {};
	template<> struct precision_cmp< double , int > : public precision_cmp_int_double {};
	// int & long double
	struct precision_cmp_int_ldouble
	{
		typedef float left;
		typedef long double right;
	};
	template<> struct precision_cmp< int , long double > : public precision_cmp_int_ldouble {};
	template<> struct precision_cmp< long double , int > : public precision_cmp_int_ldouble {};
	// float & double
	struct precision_cmp_float_double
	{
		typedef float left;
		typedef double right;
	};
	template<> struct precision_cmp< float , double > : public precision_cmp_float_double {};
	template<> struct precision_cmp< double , float > : public precision_cmp_float_double {};
	// float & long double
	struct precision_cmp_float_ldouble
	{
		typedef float left;
		typedef long double right;
	};
	template<> struct precision_cmp< float , long double > : public precision_cmp_float_ldouble {};
	template<> struct precision_cmp< long double , float > : public precision_cmp_float_ldouble {};


	// precision select
	template< typename T1 , typename T2 > struct precision_select
	{
		typedef typename precision_cmp< T1 , T2 >::right result;
	};

	// recursive op
#define	SELECT_TYPE2( OP , T1 , T2 )				typename OP< T1 , T2 >::result
#define	SELECT_TYPE3( OP , T1 , T2 , T3 )			typename OP< T1 , SELECT_TYPE2( OP , T2 , T3 ) >:: result
#define	SELECT_TYPE4( OP , T1 , T2 , T3 , T4 )		typename OP< T1 , SELECT_TYPE3( OP , T2 , T3 , T4 ) >::result
#define	SELECT_TYPE5( OP , T1 , T2 , T3 , T4 , T5 )	typename OP< T1 , SELECT_TYPE4( OP , T2 , T3 , T4 , T5 ) >::result

#undef max
#undef min
#define	max max
#define	min min

	/* comparison - extend to 6 parameters */
	template< typename _Ty > inline _Ty max( const _Ty a , const _Ty b ) { return __max( a , b ); }
	template< typename _Ty > inline _Ty min( const _Ty a , const _Ty b ) { return __min( a , b ); }

	// clamp value & return
	template< typename _Ty > static inline const _Ty& clamp( const _Ty& value , const _Ty& minimum , const _Ty& maximum ) {
		if( value < minimum ) return minimum;
		if( value > maximum ) return maximum;
		return value;
	}

	// clamp reference
	template< typename _Ty > static inline void clamp_ref( _Ty& value , const _Ty& minimum , const _Ty& maximum ) {
		if( value < minimum ) value = minimum;
		else if( value > maximum ) value = maximum;
	}

	// map unit ([0,1]) to range ([0,range])
	template< typename _Ty > static inline _Ty map_unit( const _Ty& value , const _Ty& range ) {
		return clamp( value , _Ty(0) , _Ty(1) ) * range;
	}

	/* fast sin cos */
	static inline void sincos( const float radians , float * sine , float * cosine );
	static inline void sincos( const double radians , double * sine , double * cosine );

	/* sin(x)/x, which is accurate for very small x. */
	template<typename _Ty> inline _Ty sinx_over_x( _Ty x ) {
		if( x * x < ArithmeticOp<_Ty>::epsilon() )
			return _Ty(1);
		else
			return sin(x) / x;
	}

	/* fast operations */
#define	sign( X )					( ( X > 0 ) ? 1 : ( X < 0 ) ? -1 : 0 )

	/* epsilon comparison */
	template<typename _Ty> inline bool is_zero( const _Ty& a ) {
		return ( fabs(a) < ( ArithmeticOp<_Ty>::comparison_epsilon() ) );
	}
	template<typename _Ty> inline bool is_equal( const _Ty& a , const _Ty& b ) {
		return ( ( fabs(_Ty( a - b )) ) < ( ArithmeticOp<_Ty>::comparison_epsilon() ) );
	}
	/*	Length of a value */
	template<typename _Ty> inline const _Ty& length( const _Ty& val ) {
		return val;
	}
	/*  Degree 2 Radians */
	template<typename _Ty> inline _Ty radians( const _Ty& angle ) {
		return angle * ( ArithmeticOp<_Ty>::pi_over_180() );
	}

	/*	Radians 2 Degree */
	template<typename _Ty> inline _Ty degrees( const _Ty& radian ) {
		return radian * ( ArithmeticOp<_Ty>::pi_under_180() );
	}

	/* reciprocal */
	template<typename _Ty> inline _Ty recip( const _Ty& a ) {
		return _Ty( 1.0 ) / a;
	}

	/* reciprocal of square root */
	template<typename _Ty> inline _Ty recip_sqrt( const _Ty& a ) {
		return _Ty( 1.0 ) / ( sqrt(a) );
	}

	/* round */
	template<typename _Ty> inline _Ty round( const _Ty& a ) {
		int i = (int)a;
		if( ( a - _Ty(i) ) < 0.5 )
			return _Ty(i);
		return _Ty( i + 1 );
	}

	/* gamma */
	template< typename _Ty > inline _Ty gamma( const _Ty& g , const _Ty& x ) {
		return pow( x , recip( g ) );
	}

	/* bias */
	template< typename _Ty > inline _Ty bias( const _Ty& b , const _Ty& x ) {
		return pow( x , log( b ) / log( _Ty(0.5) ) );
	}

	/* gain */
	template< typename _Ty > inline _Ty gain( const _Ty& g , const _Ty& x ) {
		if( x < 0.5f )
			return bias<_Ty>( _Ty(1) - g , _Ty(2) * x ) / _Ty(2);
		else
			return _Ty(1) - bias( _Ty(1) - g , _Ty(2) - _Ty(2) * x ) / _Ty(2);
	}

	/* step */
	template<typename _Ty> inline _Ty step( const _Ty& a , const _Ty& b ) {
		return b < a ? _Ty(0) : _Ty(1);
	}

	/* smooth step */
	template<typename _Ty> inline _Ty smoothstep( const _Ty& min , const _Ty& max , const _Ty value ) {
		if( value < min )
			return 0;
		else if( value > max )
			return 1;
		else
		{
			_Ty u = ( value - min ) / ( max - min );
			_Ty v = u * u;
			return -2.f * u * v + 3.f * v;
		}
	}

	/* lerp */
	template<typename _Ty> inline _Ty lerp( const _Ty& a , const _Ty& b , const float t ) {
		return a + b * t;
	}

	/* mix */
	template<typename _Ty> inline _Ty mix( const _Ty& a , const _Ty& b , const _Ty& alpha ) {
		return ( a * ( 1 - alpha ) ) + ( b * alpha );
	}

	/* fast blend with [0,255] value */
	inline byte blend( const uint src , const uint dst , const uint alpha )
	{
		return byte( ( ( ( src - dst ) * alpha >> 8 ) + dst ) & 0xff );
	}

	/* fast blend with [0,255] value */
	inline byte blend_alpha( const uint src , const uint dst )
	{
		//alpha = (unsigned char) ((alpha+(255-alpha)*(al/255.f)));
		return byte( ( dst + ( ( ( 255 - dst ) * src ) >> 8 ) ) & 0xff );
	}

	/* spline */
	template< typename _Ty > inline _Ty spline( const _Ty basis[16] , const _Ty value , const _Ty knot[4] );

	/* polynomial solving */
	template <class _Ty> inline int solveLinear( const _Ty& a , const _Ty& b , _Ty &x );
	template <typename _Ty> inline int solveQuadric( const _Ty& a , const _Ty& b , const _Ty& c , _Ty r[2] );
	template <typename _Ty> inline int solveCubic( const _Ty c[4] , _Ty s[3] );
	template <typename _Ty> inline int solveQuartic( const _Ty c[5] , _Ty s[4] );

	/* (2d triangle area ) * 2 */
	template<typename _Ty> inline _Ty area( _Ty x1 , _Ty y1 , _Ty x2 , _Ty y2 , _Ty x3 , _Ty y3 )
	{
		const double acx = x1 - x3;
		const double bcx = x2 - x3;
		const double acy = y1 - y3;
		const double bcy = y2 - y3;
		return _Ty(acx * bcy - acy * bcx);
	}

	/* adjoint of matrix */
	template< typename _Scalar > inline void adjointx( _Scalar *out,const _Scalar *in );

/* vector arithmetic operations */

// dot product
#define	dot2( a , b )			( ( (a)[0] * (b)[0] ) + ( (a)[1] * (b)[1] ) )
#define	dot3( a , b )			( ( (a)[0] * (b)[0] ) + ( (a)[1] * (b)[1] ) + ( (a)[2] * (b)[2] ) )

// cross product
#define	cross3( c , a , b )								\
	{													\
		(c)[0] = ( (a)[1] * (b)[2] - (a)[2] * (b)[1] );	\
		(c)[1] = ( (a)[2] * (b)[0] - (a)[0] * (b)[2] );	\
		(c)[2] = ( (a)[0] * (b)[1] - (a)[1] * (b)[0] );	\
	}

// interpolate
#define	interpolate2( c , a , b , val )			\
	{											\
		(c)[0] = (a)[0] * (1-(val)) + (b)[0] * (val);		\
		(c)[1] = (a)[1] * (1-(val)) + (b)[1] * (val);		\
	}
#define	interpolate3( c , a , b , val )			\
	{											\
		(c)[0] = (a)[0] * (1-(val)) + (b)[0] * (val);		\
		(c)[1] = (a)[1] * (1-(val)) + (b)[1] * (val);		\
		(c)[2] = (a)[2] * (1-(val)) + (b)[2] * (val);		\
	}

//	length
#define	length2( a )		( sqrt( dot2( a , a ) ) )
#define	length3( a )		( sqrt( dot3( a , a ) ) )

/* matrix related operations */

// det 2 x 2
#define det2x2( a , b , c , d ) 	 (double(a) * double(d) - double(b) * double(c))

// det 3 x 3
#define det3x3( a1 , a2 , a3 , b1 , b2 , b3 , c1 , c2 , c3 ) \
			( (a1) * det2x2( b2 , b3 , c2 , c3 ) - (b1) * det2x2( a2 , a3 , c2 , c3 ) + (c1) * det2x2( a2 , a3 , b2 , b3 ) )
// det 4 x 4
#define det4x4( m ) \
		( double(m[0]) * det3x3( m[5] , m[9] , m[13] , m[6] , m[10] , m[14] , m[7] , m[11] , m[15] ) \
		- double(m[1]) * det3x3( m[4] , m[8] , m[12] , m[6] , m[10] , m[14] , m[7] , m[11] , m[15] ) \
		+ double(m[2]) * det3x3( m[4] , m[8] , m[12] , m[5] , m[9 ] , m[13] , m[7] , m[11] , m[15] ) \
		- double(m[3]) * det3x3( m[4] , m[8] , m[12] , m[5] , m[9 ] , m[13] , m[6] , m[10] , m[14] ) )


	/* SSE optimized intrinsic */

	/* a fix of cvtsd_f64 & cvtss_f32 */
#ifndef	__ICL
	inline double _mm_cvtsd_f64( __m128d a );
	inline float _mm_cvtss_f32( __m128 a );
#endif

	/* component access */
	template< int comp > inline __m128 _mm_comp_pd( __m128d a );
	template< int comp > inline __m128 _mm_comp_ps( __m128 a );
	inline __m128d _mm_sum_pd( __m128d x );

	/* x[0] + x[1] + x[2] */
	inline __m128 _mm_sum_ps3( __m128 x );

	/* x[0] + x[1] + x[2] + x[3] */
	inline __m128 _mm_sum_ps4( __m128 x );

	/* x[0] * x[0] + y[0] * y[0] */
	inline __m128d _mm_square_pd( __m128d x );

	/* x[0] * x[0] + y[0] * y[0] + z[0] * z[0] */
	inline __m128 _mm_square_ps3( __m128 x );

	/* x[0] * x[0] + y[0] * y[0] + z[0] * z[0] + w[0] * w[0] */
	inline __m128 _mm_square_ps4( __m128 x );

	/* x[0] * x[1] + y[0] * y[1] */
	inline __m128d _mm_dot_pd( __m128d x , __m128d y );

	/* x[0] * x[1] + y[0] * y[1] + z[0] * z[1] */
	inline __m128 _mm_dot_ps3( __m128 x , __m128 y );

	/* x[0] * x[1] + y[0] * y[1] + z[0] * z[1] + w[0] * w[1] */
	inline __m128 _mm_dot_ps4( __m128 x , __m128 y );

	/* cross */
	inline __m128 _mm_cross_ps3( __m128 a , __m128 b );

	/* r = 1 / a */
	/* [2 * rcpps(x) - (x * rcpps(x) * rcpps(x))] */
	inline __m128 _mm_rcpnr_ps( __m128 a );
	inline __m128 _mm_rcpnr_ss( __m128 a );

	/* r = 1 / sqrt(a) */
	/* 0.5 * rsqrtss * (3 - x * rsqrtss(x) * rsqrtss(x)) */
	inline __m128 _mm_rsqrtnr_ps( __m128 a );
	inline __m128 _mm_rsqrtnr_ss( __m128 a );

	/* abs */
	inline __m128 _mm_abs_ps( __m128 a );
	inline __m128d _mm_abs_pd( __m128d a );

	/* negative */
	inline __m128 _mm_neg_ps( __m128 a );
	inline __m128d _mm_neg_pd( __m128d a );

	/* comparison with epsilon */
	inline bool _mm_equal_ps4( __m128 a , __m128 b );
	inline bool _mm_equal_ps3( __m128 a , __m128 b );
	inline bool _mm_equal_pd( __m128d a , __m128d b );

	/* reflect */
	inline __m128 _mm_refl_ps( __m128 n , __m128 i );

	/* refract */
	inline __m128 _mm_refr_ps( __m128 n , __m128 i , __m128 eta );

	/* normalize */
	inline __m128 _mm_normalize_ps( __m128 a );

	/* clamp value */
	inline __m128 _mm_locut_ps( __m128 val , __m128 bound );
	inline __m128 _mm_hicut_ps( __m128 val , __m128 bound );
	inline __m128 _mm_clamp_ps( __m128 val , __m128 min , __m128 max );

	/* pack vector */
	inline __m128 _mm_pack_ps3( __m128 a , __m128 b , __m128 c );
	inline __m128 _mm_pack_ps3( float a , float b , float c );
	inline __m128 _mm_pack_ps4( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_pack_ps4( float a , float b , float c , float d );

	/* select op */
	inline __m128 _mm_select_eq_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_lt_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_le_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_ge_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_neq_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_nlt_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_nle_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_ngt_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );
	inline __m128 _mm_select_nge_ps( const __m128& a , const __m128& b , const __m128& c , const __m128& d );

	/* comparison boolean */
	inline __m128i _mm_bcmp_eq_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_lt_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_le_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_gt_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_ge_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_neq_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_nlt_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_nle_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_ngt_ps( __m128 a , __m128 b );
	inline __m128i _mm_bcmp_nge_ps( __m128 a , __m128 b );

	/* select op */
	inline __m128d _mm_select_eq_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_lt_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_le_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_ge_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_neq_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_nlt_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_nle_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_ngt_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );
	inline __m128d _mm_select_nge_ps( const __m128d& a , const __m128d& b , const __m128d& c , const __m128d& d );

	/* comparison boolean */
	inline __m128i _mm_bcmp_eq_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_lt_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_le_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_gt_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_ge_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_neq_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_nlt_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_nle_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_ngt_ps( __m128d a , __m128d b );
	inline __m128i _mm_bcmp_nge_ps( __m128d a , __m128d b );

	/* sign op */
	inline __m128 _mm_sign_ps( __m128 a );

	/* fast reciprocal */
	template<>inline float recip<float>( const float& a );
	template<>inline double recip<double>( const double& a );

	/* fast reciprocal square root */
	template<>inline float recip_sqrt<float>( const float& a );
	template<>inline double recip_sqrt<double>( const double& a );

	static inline __m128 _mm_floor_ps( __m128 m );
	static inline __m128 _mm_ceil_ps( __m128 m );
	static inline __m128 _mm_round_ps( __m128 m );
	static inline __m128 _mm_radians_ps( __m128 m );
	static inline __m128 _mm_degrees_ps( __m128 m );
	static inline __m128 _mm_mod_ps( __m128 m , __m128 f );

}  // namespace xeres

#include "lib/math/math_ext.inl"

using xeres::min;
using xeres::max;
